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ABSTRACT 

We present a new computational scheme aimed at reducing the complexity of the 
chemical networks in astrophysical models, one which is shown to markedly improve 
their computational efficiency. It contains a flux-reduction scheme that permits to deal 
with both large and small systems. This procedure is shown to yield a large speed-up 
of the corresponding numerical codes and provides good accord with the full network 
results. We analyse and discuss two examples involving chemistry networks of the 
interstellar medium and show that the results from the present reduction technique 
reproduce very well the results from fuller calculations. 

Key words: astrochemistry - ISM: evolution, molecules - methods: numerical. 



1 INTRODUCTION 

Chemistry is one of the most impo rtant ingredients in 
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cesses, usually represented through a network, consist of a 
large ensemble of reactions that tightly connect a given set 
of chemical species. From a mathematical point of view this 
network is represented by a system of ordinary differential 
equations (ODE) that often constitutes a system of stiff, 
coupled equation. Computationally speaking, it means that 
a lot of the cpu-time is spent to solve the chemical network 
evolution instead of solving, for instance, the hydrodynam- 
ical equations. 

To deal with this problem there are various tech- 
niques which can be applied. The most used is probably 
the direct reduction approach that consists in pre-selecting 
the most important species and then their most impor- 
tant reactions for a ce rtain astrophysical environment (e.g. 
iNelson fc LangeJI 19991 ). This method is trivial when the net- 
work is small, or when there are small regions of the net- 
work that are easy to visualize and to isolate, but when the 
complexity grows this approach could result in large errors, 
since the network reduction is based on ad hoc chemical or 
physical considerations applied only once before solving the 
ODEs. 

Another widely used approach is to reduce the num - 
ber of non-equilibriurrQ species (e.g. iGlover et al.] l2010h . 
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Under this assumption the differential equations are ex- 
plicitly solved only for some selected species while the re- 
maining species are taken to reach instantaneously their 
equilibrium. Even this reduction technique is strongly 
problem-dependent and this approximation could lead to 
large uncertainties in the final abundances. 

The aim of the present work is to suggest instead a 
method that increases the efficiency of the system of differ- 
ential equations which describes the chemical network, but 
without introducing any a priori assumption, while provid- 
ing an approximate but still accurate solution that matches 
that from the full network calculation. 

The general method was developed by iTuppeil (|2002l ) 
(hereafter TUP02) while we introduce modifications to fit 
our class of astrophysical scenarios. We shall show that the 
TUP02 scheme with our improvements yields accurate re- 
sults while achieving large computational speed-ups. We 
shall a pply the method t o a small chemical network (sim- 
ilar to lGlover et al.ll201(]| ) and to a larger one (similar to 
IWakelam fc Herbstll20ol ). We first introduce the mathemat- 
ical problem in Section [2l then we describe in Section [3] the 
main features of the present scheme and show our modifica- 
tions and upgrades to the original method. Finally, Section 
[4] discusses the chosen astrophysical examples. 



2 EARLIER NETWORK REDUCING 
METHODS: AN OUTLINE 

A chemical network with A'^ species is represented by a sys 
tem of A'^ ordinary differential equations as 
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where i G K^, rii is the abundanc^l of the ith specie^ and 
kip is the reaction rate that represents the probability for 
the reaction between the ith and the j'th species to occur at 
a chosen temperature. The right-hand side (RHS) of Eq.([l]) 
is divided into two parts: the first is a sum that describes 
the reactions which destroy the ith species while the second 
part is a sum over the reactions that produce the ith species. 
Here i, p, q, and r indicate each a specific species. 
Eq.(IT} can be now written as 

M 

hi =^SijFj{n,k) , (2) 

where Sij is an element of the stoichiometric matrix, 
Fj = kimniUm is a fiux, and M is the total number of reac- 
tions. The elements of the stoichiometric matrix are integers 
and can be Sij — ±1 (depending on the direction of the reac- 
tion) and Sij — when the reaction does not occur. The fiux 
is the amount of mass or number density that moves from 
a species to another in a unit of time. Equation Q allows 
to viewing the network as a dtrected graph, where each node 
(or vertex) represents the number density of each species, 
and each link (or edge) represents a flux. The direction (and 
the existence) of an edge is determined by the correspond- 
ing element of the stoichiometric matrix s. Note that here i 
indicates a species, while j now represents a flux. 

This system can be easily solved with DVODE 
ijHindmarsh' 1983h using a backward differentiation formulas 
(BDF) integration scheme. The computational cost of solv- 
ing the ODE system is given by the evaluation of the RHS 
of the Eq.((2]l and in some cases also by the large number 
of equations. Hence we need to find a way of reducing the 
complexity of the system. There are two possible approaches 
to this end: one is based on reducing the number of differ- 
ential equations and another keeps the number of equations 
constant but reduces the number of elements in the RHS of 
each equation. These methods are not mutually exclusive. 

The first technique to reduce the system consists in re- 
duci ng the number of species throug h some lumping method 
(e.g. lOkino fc Mavrovouniotiall998l ) like the singular value 
decomposition. In this case the full system - see Eq.Q - 
becomes a system oi K < N equations 

^=^J.,i^(n,fc), (3) 

where i £ R^. The hat indicates the same quantities of 
Eq.((2} but transformed by a unitary N x K matrix U that 

permits R^ and analogously its transposed U' lets 

R'^ <^ R^'. This method allows to solve the ODE system 
in a different space where h £ R^ are linear combinations 
of the various n G R and one chooses the transformation 
U in order to obtain a simpler system. Unfortunately the 
main drawback of species-based reduction is that building 
the transformation matrix U may not be trivial and also 
the computational cost to build U could easily result in a 
computational overhead. 

^ Abundances are always considered as number densities unless 
otherwise indicated. 

^ A species can be both an element or a molecule, e.g. H2 and Si 
are both species. 



In the present work we therefore follow a different ap- 
proach that consists in reducing the ffuxes between the 
species while keeping the same number of differential equa- 
tions. In other words, we look for a system with W < M 
represented by the following ODE system 

w 

hi = '^SijFj{n,k) , (4) 

j=i 

where i G R^, but in this case the sum runs over fewer el- 
ements than those of the system of Eq.Q. The advantage 
of this method is that the greater part of the computational 
cost is given by the RHS of Eq.([2]| and we reduce the num- 
ber of reactions instead of reducing the number of species. In 
particular, we delete the reactions that are slow or, in other 
words, the reactions that have small fiuxes. This method 
has some practical advantages compared to the one which 
is species-based: (i) the mass conservation is always guar- 
anteed because we are removing whole reactions, and (ii) 
the selection of small fluxes is more efficient than building a 
unitary transformation, thus decreasing the computational 
overhead. On the other hand, the drawback of this method 
is that it requires a robust selection criterion to determine 
what are the important reactions. 

A combined speci es-based and reacti on-based methods 
has been sug gested bv lWiebe etHI (|2003l ) and we focus our 
attention on their reaction-based method because it is more 
interesting for the present discussion. They search for the 
production and destruction reactions that are most impor- 
tant for the evolution of an important species, and then 
they assign their own relative significance. The most impor- 
tant species are determined by an iterative algorithm that 
computes an importance weight for each species to deter- 
mine what reactions are important or not. Their interesting 
method led to good results but unfortunately it seems to 
be problem-dependent and, moreover, it is affected by un- 
expected changes of the physical conditions of the evolving 
model. The latter is the most troublesome d rawback. 

The method proposed bv lTupperl l|2002h belongs to the 
class of the reaction-based methods: the main assumption 
is to consider only the fiuxes that are larger than a given 
default value, namely e, in order to select the fiuxes only 
if £ < \Fi\, this procedure being repeated at each macro- 
step, instead of doing it once for all at the beginning of the 
simulation. It is called macro-step to underline the fact that 
is larger than the usual integration steps, so that a macro- 
step includes o ne or more integrations steps. This is based on 
the approach of lPetzold fc Zhul l|l999l ) and already described 
by Eq.Q, although TUP02 introduces some differences, the 
most important of which consist of using different reduced 
models over short periods of time instead of a single reduced 
system over the entire integration interval. This approach 
turns out to be more effective although problems may arise 
since: (i) one needs to determine the size of a macro-step 
and (ii) the method requires an error control to check the 
validity of the macro-step integration. 

To solve the first problem TUP02 considers the infiu- 
ence of the omitted reactions to predict for how long the 
model will remain valid. In particular one needs to know 
the infiuence of the omitted reactions over the next macro- 
step or, even better, when such reactions will become large 
enough to make the chosen macro-step not valid. TUP02 
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uses finite differences to compute tfie rate of change of a sin- 
gle reaction obtaining Gj ~ dFj{n) /dt. Then he considers 
that the reaction most likely to cause a macro-step rejection 
is the one with the largest Gj, this assumpti on permitting 
to find the size of a macro-step (|Tuppeill200'3 . Section 4.3). 

The second problem depends on the omitted reactions: 
if they become too large (larger than a default tolerance) 
during a single macro-step, the integration over that time 
period is not valid. To ensure the effectiveness of the reduc- 
tion TUP02 suggests to compare the omitted reactions at 
the beginning and at the end of the macro-step. If there is 
any difference between these two check-points the macro- 
step is rejected and then the entire integration is repeated 
including the omitted reactions. This procedure is iterated 
until the difference disappears. Our present scheme, how- 
ever, does not use any rejection technique because we use a 
more accurate macro-step, as will be discussed in Section [3l 
so that for our models we found that a rejection scheme is 
not necessar y. 

iTuppeil l|2002l ) then applies the method to a large as- 
trophysical network and to a small combustion network. For 
the first case the speed-up is very large, and the differences 
of results between the full model and the reduced one are 
negligible. In the combustion case the full system is well re- 
produced, but it becomes computationally inefhcient since 
a large overhead occurs, so that the reduction method is 
slower than the straightforward integration. It shows that a 
massively interconnected system could be difficult to reduce, 
in particular when all the fluxes have the same absolute val- 
ues. On the other hand, in the first example, the algorithm 
easily finds a subset of reactions that is dominant compared 
to the others. 



3 THE PRESENT REDUCTION SCHEME 

Our a pproach is similar to the one proposed by iTuppeij 
although we introduce the following important mod- 
ifications: 

1. Adaptive tolerance: we define the tolerance e rel- 
ative to the largest flux as £ = |-Fmax|. This relative toler- 
ance changes during the evolution of the system, depending 
on the magnitude of the largest flux, the one that domi- 
nates the model. We then select the fluxes that obey to the 
following rule 



\F,\>e = (;\Fn- 



(5) 



where depends on the desired accuracy (the best value 
for will be discussed later). 

2. Ma cro-ste p leng th: another modification to the 
scheme of iTupped (|2002l ) is the choice of the macro- 
step length. We first consider the aforementioned equation 
Gj « dFj{n)/dt. This equation can be viewed as the veloc- 
ity of the jth flux: in other words, it shows how much the 
flux of a given reaction changes during a time interval. We 
use this quantity to predict when a reaction will cross the 
value e. Figure [T] represents the amount of different fluxes, 
with the line labelled e indicating the tolerance value. The 
fluxes in the hatched area above that line obey the condi- 
tion while the fluxes below the e-line are not considered 
in the reduced system. Each horizontal solid line is the flux 
value at a given time t, and the corresponding dashed line 



B 



D 



A 



IE 



Figure 1. The evolution of fluxes in a macro-step. See text for 
further details. 



is the value of the same flux at time t + At, as indicated by 
the arrow which represents the variation of flux over At; the 
larger the Gj the longer the arrow. We show five different 
cases: (A) a reaction important at time t becomes uninter- 
esting after At, since it crosses the e-line to reach the darker 
area, (B) a reaction becomes important after At, (C) a reac- 
tion for which the contribution remains important, and (D, 
E) two reactions which remain unimportant even if one of 
them increases its importance in a macro-step. To estimate 
the length of a macro-step we are interested in the reactions 
that belong to the classes (A) and (B), since their fiuxes 
change their significance before the end of a macro-step. We 
calculate the macro-step length as 
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At 



(6) 



where {F^ — e) represents the jth fiux distance from the 
threshold e, and {Fj — i?^*+'^*) is the fiux variation over the 
previous macro-step of length At. A small value is intro- 
duced to avoid divisions by zero. The constant factor (j)^ 1 
is the time-tolerance, a parameter to avoid excessively short 
macro-steps caused by the reactions that oscillate around 
the threshold e. Note, however, that if cf> is chosen to be too 
large it can affect the accuracy of the method (in this paper 
we use (f) = 10). 

It is worth noticing that the direction of the flux velocity 
is important, and then fast growing reactions that are above 
the threshold should be ignored when one determines the 
length of a macro-step. This consideration is also true for 
the fast decreasing reactions that are smaller than e. To take 
into account this feature and to avoid shorter, unnecessary 
macro-step lengths in Eq.(|6}, one should write it as 



Atno 



pt _ pt + At 
3 3 



At 



(7) 



although for programming purposes we shall still use Eq.®. 
The reader must be aware that this macro-step calculation 
is based on a linear growth hypothesis (see Gj) and it could 
be use d together with a rejection scheme, as done bv lTuppeil 
(j200^) , if the system presents a strong non-linear behaviour. 
Again it should be mentioned that for our modelling we 
found it to be unnecessary. 
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Figure 2. Small model calculations. Time evolution of C+ (top) 
and OH (bottom) number densities. Open circles indicate the full 
model evolutions, while lines are for different C, in the model. All 
lines are clearly overlapping. 



3. Sorting: We further improved on the TUP02 
method by introducing a flux-sorting device at each macro- 
step. The system of Eq.© is built with a loop over all the 
M reactions, and each reaction is evaluated with the con- 
dition ([S]), hence we must evaluate this condition M times: 
this operation has a non- negligible computational cost. If 
the fluxes are sorted out, then we can break the loop after 
W iterations - see Eq.@. This strategy is not effective if the 
macro-steps become very small, because the time to sort the 
fluxes becomes larger than the time gained in the early loop- 
breaking. In the astrophysical tests employed in the present 
paper we never incurred in a computational overhead, even 
if we used a simple bubble sorting. One should note here 
that the first sorting is the most expensive while the others 
are very efficient, because they find an almost-sorted array. 
This reduces the inefficiency of the bubble sorting proce- 
dure (|Knuthll 19971 ). To increase the sorting efficiency we in- 
troduce a boolean array where each element is false or true 
depending on condition ([Sjl . We sort then the boolean array 
instead of the array that contains the fluxes because we are 
only interested in evaluating the condition ((5)l rather than 
generating sorted fluxes. 

The overall efficiency of our approach is determined by 
the time spent to prepare a macro-step Tm compared to the 
solver's integration time over a macro-step. For a macro-step 
Atmacro we must have 



where rdo is the time to calculate a single do-cycle-step, 
Afodc is an integration step of the ODE solver. The last ra- 
tio determines the number of ODE solver's calls in a macro- 
step to the function that computes the RHS of Eq.((2|. We 
roughly assume that without sorting the cost of creating a 
macro-step is given by Tm = 3M (namely M to find the max- 
imum, M to create the boolean array, and M to compute 
the length of the macro-step, where M is defined in Eq.((2| 
in the previous Section) when introducing the bubble sort- 
ing, the worst possible case is Tm ~ 3M + . Analogously, 
the reduced do-loop needs Tdo = Wrdn + Mth steps to cal- 
culate the RHS of our differential equations, where W is 
defined in Eq.Q, rdn is the cost of evaluating An, and m is 
the cost of evaluating the condition ((5]). When we introduce 
sorting we have th — 0. The efficiency is then guaranteed if 
C > 3A//(Wrdn + Mrif) andffC > {3M + M'^)/(WTin) with 
and without sorting, respectively. It is worth to mention at 
this point that using sorting for our examples enabled us to 

obtain Taorted ~ 0. 75 Tunsorted . 

We illustrate our scheme by briefly outlining a pseudo- 
code (see Algorithm [T| , where n is an array containing the 
abundances of the species, t is the simulation time, tmacro is 
the length of a macro-step, ttarg is the target time to solve 
the system, t^nd is the length of the simulation, e is the flux 
limit, is deflned in Eq.([5ll, F is the array of the fluxes, and 
u is the aforementioned boolean array. 

As an implementation of our scheme, this pseudo-code 
will be applied to two different astrophysical scenarios in the 
following Section. 



4 RESULTS AND DISCUSSION 
4.1 The small network example 

We now apply our scheme to a one-zone che mical network 
using a set of reactions similar to the one in I Glover et al.l 
involving 29 species and 170 reactions (see Tab lAll 
for the complete list). The chosen subset of reactions has 
only an illustrative purpose since we are not interested in 
reproducing the physical behaviour of a real astrophysical 
environment, but we want only to test the reduction method 
proposed. Hence, this set of reaction is intended only to 
mimic an astrophysical behaviour. These considerations are 
also true for the large net work example. 

We also note that in iGlover et al.l (|2010l ) thirteen of 
the original species are considered to be in instantaneous 
chemical equilibrium, while the remaining nineteen species 
follow the full non-equilibrium evolution. In our model none 
of the 29 species are assumed to be in instantaneous equi- 
librium so that we track the non-equilibrium evolution for 
all of them. Our initial conditions are n = 10~^° cm~^ for 
all the species except tih = lO^cm""^, no = 3.16 x 10~*nH, 
and nc = 1.41 x lO'^^nn. We also set A„ = 10, and T = 10^ 
K, then we let the system evolve until tend = 10^ yr. For our 
test we did not use any cooling or heating, hence dT/dt = 0. 

The program is serial and it is written in FOR- 
TRAN 90, compiled with Intel® Fortran Compiler 12. iQ on 
an Intel® Xeon® E5430. The solver used is the DVODE 



Tm < Tdo- 



Atn 



Atod 



= Tdo ■ C , 



(8) 



* Optimization flag used: -03, 
details contact the authors. 



unroll, -ip, and -ipo. For further 
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Algorithm 1 - Pseudo-code of the scheme proposed in this work. See text for further details. 



n no 
t ^ 

imacro ^ 0.01 
^targ ^ ^macro 

repeat 

max(F*) 
u -f- evaluate (e, 
[m, -F] sort(M) 

tmacro get_macrostep_length( AF, At) 

ttarg max(f + t macro 5 ^cnd) 
odesolve(f, ttarg, n, F, u) 
until (t < tend) 



/ /initialize n 
I /initialize t 

1 1 set a default value for the first macro-step 

/ /set the solver target time 

/ /loop over the macro-steps 

/ / compute flux limit 

/ /make boolean array 

/ / sort u and F using u 

1 1 compute the length of the next macro-step 

/ /set next target time 

/ /solve the reduced system for [t, ttarg] 

/ / exit when the end is reached 



l|Hindmarshlll983h with absolute tolerance 10 *° and rel- 
ative tolerance 10~^. 

We have made flve simulations for flve different values 
of C^. Our aim is to compare the time evolution of the various 
species provided by the full model with the evolution from 
the reduced ones. The larger fraction of the species behaves 
as the examples shown in Fig[2] and Fig[3] where the full 
model is well reproduced by all the reduced models we em- 
ploy. Some of the species evolve as in Fig[4l in this case the 
full model starts to be well reproduced only when the impor- 
tance of these species crosses a given value that depends on 
C,. We note that different reduced models eventually catch- 
up with the full model at different times, depending on the 
selected C,. It is important to remark that only the species 
with an overall small number density are affected by this 
behaviour. 

The computational speed-ups are indicated in Tab[T] 
We also see that lower values of C, imply slower integration 
times because, when the threshold is small, there are more 
reactions that affect the construction of the RHS of Eq.Q. 
The number of reactions used at a given time is shown in the 
upper panel of Fig[S] For all the reduced models the number 
of reactions grows with time, because more reactions become 
important in the network when new species are formed. The 
number of reactions used at any given time is also propor- 
tional to since more reactions cross the threshold e, also 
proportional to ^. This is consistent with noting that, when 
C = 0, the reduced model corresponds to the full model and 
- analogously - when ^ — >■ the reduced models converge 
to the full one. Hence, when C, becomes very small the full 
model is almost exactly reproduced but a computational 
overhead would occur. Using Eq.([8]) it is possible to deter- 
mine the value ("o that maximizes the accuracy without over- 
heads, so when C < Co the reduction method is efficient and 
the choice of depends on the desired accuracy /efficiency 
trade-off. These considerations suggest that the user must 
choose the value of before doing the calculation, depend- 
ing on the necessary accuracy and the efficiency. However, 
the tests proposed show that = 10~^ represents a good 
compromise. 

Note that (Tab[T]) we obtain better timing when we in- 
troduce sorting, as expected, except for the large network 
when C, — 10~^. This unexpected result could be caused 
by a cache-locality effect which is architecture-dependent as 
proved by tests made on different machines (not shown here 
because the details are not relevant). It is also important to 
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Figure 3. Small model calculations. Time evolution of CO (top) 
and H2 (bottom) number densities. Open circles indicate the full 
model evolutions, while the lines are for different f . All lines are 
clearly overlapping. 



note that without using the sorting procedure, for the small 
network example with Q = 10"^'^, an overhead occurs. Con- 
versely, introducing a sorting procedure allows not only to 
avoid such overhead, but also to obtain a large speed-up. 

Using our present reduction method we can also check 
the importance of the different reactions in a given network 
by simply measuring the time that a chosen reactions is 
used by a reduced network. We compare these times with 
the total integration time tend employed for the model with 
C = 10"" and we show the results for this comparison in 
the last column of Tab lAll Note that the first 113 reactions 
play a role in the model for more than 50% of the total time. 
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Figure 4. Small model calculations. Time evolution of H+ (top) 
and CH^ (bottom) number densities. Open circles indicate the 
full model evolutions, while lines are for different C,. 
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Figure 5. The time evolution of the reactions used in the reduced 
model with different C, values. The horizontal line indicates the 
full model evolutions (constant reaction number), while varying 
lines are for different ^. Top panel: small model, bottom panel: 
large model. Note the logarithmic scale. 



Table 1. Speed-ups for different C, and for small and large network 
examples, with and without sorting. Normalized to the full model. 
See the discussion in the text. 





small 


large 




logiO 


no sort 


sort 


no sort 


sort 


full 


1.00 


1.00 


1.00 


1.00 


-7 


0.71 


0.52 


0.07 


0.05 


-8 


0.78 


0.58 


0.08 


0.06 


-9 


0.81 


0.61 


0.08 


0.08 


-10 


0.85 


0.65 


0.16 


0.11 


-11 


0.88 


0.68 


0.20 


0.14 


-12 


0.99 


0.71 






-13 


1.04 


0.75 







while the remaining reactions have a considerably smaller 
role over time. Table lAll is referred to the small network 
only. 

4.2 The large network example 

In this subsection we apply the present s cheme to a chemical 
netwo rk similar to the one employed bv lWakelam fc HerbstI 
l|2008l ) - see the considerations made for the small example 
at the beginning of Sect l4.1l This network is less connected 
compared to the previous one, even if there are more species 
(451) and more reactions (4399). For such a system we ex- 
pect to obtain a larger speed-up, because, due to its lower 
connectivity and then the possibility of deleting slower reac- 
tions, a set of almost-isolated subsystems could be identified. 



One should note that slow reactions are represented in our 
scheme by small fluxes, hence they can be ignored if they 
lie under the e thre shold. We follow the initi al condition 
of the EA2 model in IWakelam fc HerbstI l|2008h (except for 
rigi-i- = 2.56 X 10~*cm~^), tihs = 10"* cm~^, while the re- 
maining species are set to 10"^" cm"'', and Uc — "Y^rnon- 
We let the system evolve until t = tend = 10'' yr, with 
CcR = 1-3 X 10"^'' s-\ T = 10 K, and AT/At = 0. The 
model is dust-free. 

The results are similar to the ones found in the previ- 
ous calculation. The time evolution of the various species can 
be divided into two classes: well-reproduced evolutions (as 
in Fig|6]), and evolutions with a catch-up behaviour (FigO. 
Note that in this latter case the model with = lO"'^' almost 
matches the full model behaviour. Here, as in the small net- 
work case, the plots show a (^-depending shape for the same 
reasons previously described. 

In this case we obtain a very large speed-up, since the 
RHS of Eq[5]is massively reduced (see FigO bottom panel). 
The number of reactions (that depends on C,) grows in the 
first year and then decreases to lie on a long plateau. There 
is also more noise compared to the small model, because the 
length of the macro-steps is considerably smaller, as we see 
when the lines are around 1 yr where the system reaches its 
largest complexity and, thus, its smallest reduction. 

To provide a table analogue to T a blAll an additional 
calculation for the IWakelam fc HerbstI (|2008h network has 
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Figure 6. Large model calculations. Time evolutions of C (top) 
and Si (bottom) number densities. Open circles indicate the full 
model evolutions while lines are for different Lines are clearly 
overlapping. 
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Figure 7. Large model calculations. Time evolutions of C3S (top) 
and S"'" (bottom) number densities. Open circles indicate the full 
model evolutions while lines are for different 



been carried out. The table is available on request from the 
authors. 



5 CONCLUSIONS 

We have proposed a new network-redu cinR method ob tained 
as an extension of that described by iTuppeij l|2002l ) . This 
reaction-based reduction approach allows one to increase the 
computational efficiency of solving the ODE system repre- 
senting the chemical network. Our main modifications in- 
volve: (i) the introduction of an adaptive tolerance, (ii) a 
new definition of the macro-step length, as in Eq.®, and 
(iii) a sorting technique to obtain larger speed-ups. We have 
applied our metho d to the chemistry o f the ISM for a small 
network (si milar to Glover et al.l 2O10l ) , and for a larger one 
(similar to IWakelam fc HerbstI l20o5) . Our results show a 
very good agreement between the time-evolution of the full 
model and that of the reduced one, providing good speed- 
ups in the first case (about 50%) and a very good one for 
the larger network (approximately 90%). The accuracy and 
the speed-ups are coupled and their trade-off depends on the 
choice for the factor 

The next step on further improvements could be that of 
introducing an intermediate check during the macro-step, in 
order to control the validity of the linear-growth assumption 
in the definition of Gj. This check allows one to avoid the 
errors that arise from a non-linear behaviour, and to verify 



the consistency of the macro-step length before the end of 
the step itself. This method could be considered an alterna- 
tive to the macro-step rejection procedure in TUP02, with 
our scheme recovering the cpu-time which would be spent 
in the rejected macro-step. Our modifications should there- 
fore provide a computational gain when strong non-linear 
behaviour occurs. 

It is worth noticing here that our approach, when com- 
pared to the other reducing techniques proposed over the 
years as already discussed in Section (2] is based on a robust 
mathematical background that ensures its applicability to a 
wide range of astrochemical problems and turns out to be 
very little problem-dependent. 

This encouraging result suggests to apply this method 
to more complicated problems such as the three-dimensional 
hydrodynamical simulations of the system's evolution in the 
physical space within the ISM, where solving the ODEs 
involving the chemistry network along that simulation be- 
comes even more cpu-demanding. 
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APPENDIX A: LIST OF CHEMICAL 
REACTIONS FOR THE SMALL NETWORK 
EXAMPLE 

The list of chemical reactions is reported in Tab. lAll 

This paper has been typeset from a 1]eX/ X^T^X file prepared 
by the author. 
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Table Al. Percentage of the total time that a given reaction partici pates in the reduced network with f = 10 The values of the 
reaction rates k chosen for the small example are from iGlover et al.. i2010i ) except those involving carbon and oxygen which are from 
OSU database {osu_01_2007). This table refers to the small network example only. 
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